Heat Conduction in two-dimensional harmonic crystal with disorder 
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We study the problem of heat conduction in a mass-disordered two-dimensional harmonic crystal. 
Using two different stochastic heat baths, we perform simulations to determine the system size (L) 
dependence of the heat current (J). For white noise heat baths we find that J ~ 1/L a with a ~ 0.59 
while correlated noise heat baths gives a ~ 0.51. A special case with correlated disorder is studied 
analytically and gives a = 3/2 which agrees also with results from exact numerics. 
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INTRODUCTION 



The problem of proving or verifying Fourier's law, 
J = -kVT where k is the thermal conductivity, in any 
system evolving through Newtonian dynamics has been 
a challenge for theorists 0, 0] • So far most studies have 
been restricted to one-dimensional systems for the sim- 
ple reason that they are easier to study through simu- 
lations and through whatever analytic methods that are 
available. Also the hope is that such studies on one- 
dimensional systems provides insights which will be use- 
ful when one confronts the more difficult (and experi- 
mentally more relevant) problem of higher dimensional 
systems. For one dimensional systems, some of the most 
interesting results that have been obtained are: (1) For 
momentum conserving nonlinear systems, the heat cur- 
rent J decreases with system size L as J ~ 1/L a where 
a < 1 3J. Thus Fourier's law (which predicts a = 1) is 
not valid. The exponent a is expected to be universal 
but its exact value is still not known. A renormalization 
group analysis of the hydrodynamic equations 0] predicts 
a = 2/3 while mode-coupling theory [5j gives a = 3/5. 
The results from simulations are not able to convincingly 
decide between either of these. (2) For disordered har- 
monic systems we get J ~ l/L a again but the exponent 
a depends on the properties of the heat baths |y, 13, 111 • 

In dimensions higher than one, there are few detailed 
studies and it is fair to say that it is totally unclear as 
to whether or not Fourier's law will hold and, if not, 
then what the value of the exponent a is. For nonlin- 
ear systems which are expected to show local thermal 
equilibrium, both the hydrodynamic approach and mode 
coupling theories predict a logarithmic divergence of the 
conductivity in two dimensions. There have been simu- 
lations by Lippi and Livi who find a logarithmic di- 
vergence but simulations on larger-size systems by Grass- 
berger et al seem to obtain a power law divergence. A 
disordered harmonic model in 2D was studied in simula- 
tions by Yang [12j who claimed that beyond some critical 
disorder one gets Fourier's law, i.e a = 1. It is doubtful 
if this claim is correct. The data in the paper seems to 



indicate J ~ 1/L 2 which is not Fourier's law. Besides, 
these simulations were done with Nose-Hoover heat baths 
and it is known that these can be problematic when ap- 
plied to harmonic systems Simulations by Hu et al 
[14| on the same model but with stochastic heat baths 
do not find a Fourier behaviour. Finally an older study 
by Poetzsch et al looked at heat conduction in a 2D 
system with both disorder and nonlinearity and they give 
some evidence for Fourier behaviour. 

In this paper we consider heat conduction in a 2D dis- 
ordered harmonic system. Let us first try to see what 
one should expect theoretically. We expect localization 
phenomena (for phonons) to play an important role. A 
renormalization group calculation [16j predicts: in ID 
all modes with > 1/L 1 / 2 are localized; in 2D all modes 
with > [log(L)]" 1 / 2 are localized; in 3D there is a finite 
band of frequencies of extended states. This is similar to 
results for electron localization with the important dif- 
ference that here the 0-^0 modes are extended even in 
ID and 2D. Also for the case of electrons, only electrons 
near the Fermi-level contribute significantly to transport 
while in heat transport all phonons contribute. From the 
localization results we expect that in 3D the current in 
a disordered harmonic system should be independent of 
system size (a = 0). In one and two dimensions it is the 
small number of low- frequency phonons (0 < C ) which 
dominate transport properties. The fact that C — » 
with increasing L immediately implies that a > 0. In ID 
it has been shown jy| that the exact value of a depends 
on the low frequency spectral properties of the bath. A 
similar calculation is not available in the 2L> case and we 
address this specific question. 

Here we present results from a detailed simulational 
study to determine the exponent a for a mass-disordered 
harmonic system. Two different kinds of stochastic baths 
are considered, one with white noise and the other with 
correlated noise. We also study a special case where the 
disorder is correlated. 

Definition of model: We consider heat conduction 
in a two-dimensional mass-disordered harmonic crystal 



described by the Hamiltonian 
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where {x^ , pij , m.y } denote the position (particle dis- 
placements about equilibrium positions), momentum and 
mass of a particle at the site (i, j). We set the masses of 
exactly half the particles to one and the remaining to two 
and make all configurations equally probable. Heat con- 
duction takes place in the x-direction and we assume that 
the ends of the system are fixed by the boundary condi- 
tions xqj — = XL x+ ij. We will assume periodic bound- 
ary conditions along the y-direction so that Xij + L y — Xij. 
The heat baths arc modeled through Langevin equations 
and thus we get the following equations of motion: 



mijXij = —Axij + x 2 j + xij-i + Xij+i + hf 



mLvjXL^j = ~4%L x j + X L;c -lj + X Lx j-l + X Lx j+l 
mijXij 4Xjj ~F Xi— \j ~F Xi-^-ij ~{~ Xij — \ ~{~ Xij^-i (1) 



(for 1 < i < L x and 1 < j < L y ), where and h^ 
denote the forces from the heat baths. We will consider 
two different models for the heat baths: 

(I) Gaussian white noise source. Thus — —"fiij + 



rjf; hf = -yx Lm j + % 
properties (77^} 



j, where the noise terms have the 



t')S ]r , {rif{t)r,f{t')) = 2T Rl S(t - t')5 jr . 

(II) Gaussian exponentially correlated source. In this 
case the bath forces have the forms: 



dt'^{t-t')xi j {t')+r l 



dt'^{t-t')x Lij {t')+i 1 



(2) 



with (vf{t)rjf,(t')) = T L j(t-t')S n ,, (ri?{t)T$V)) = 
TiClit — t')5jji and ^(t) = e~ 7 *. A simple way of imple- 
menting correlated baths in the simulations is by intro- 
ducing new dynamical variables , y^ for the bath and 
setting hj — yj, h^ — y^. These satisfy the equations 
of motion yf = — TS/^ ™ %ij + Vi> e ^ c - I R t ne long time 



limit it can be easily seen that the solutions yf(t) and 
yf(t) have the required properties of correlated baths. 

We now discuss the results from simulations of the two 
different bath models. 

Simulations with white noise: Equilibration times 
in simulations of disordered harmonic lattices can be very 
long and this can sometimes lead to wrong conclusions 
[see for e.g |17j ] . To avoid such problems we first compare 
our simulation results with exact numerical results on 
steady state properties of small systems. We now briefly 
describe the numerical technique. 




FIG. 1: Temperature at all the sites of a 8 x 8 fully disordered 
lattice, from simulations and from the exact solution. Inset 
shows the disorder-averaged temperature profiles for different 
system sizes and seems to approach a linear form. 



With TV = L x L y let us define the new variables 
{<7i,<?2, -,«2Jv} = {xii,x 12 ...,x LscLy ,p 11 ,p 12 ...,p L:rLy }. 
Then Eq. ^ can be rewritten in the form: 
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(3) 



where the vector £ has all elements zero except £,N+j = 
Vf; &N-L y +j = nf (for 1 < j < L y ) and the 2N x 2N 
matrix a is given by: 

-M- 1 

$ r 

where the N x N matrices M, 4>, T can be labeled by the 
double indices (i, j) and are given by 



M, 



— 5a/(5jji-i + Sjji + i) — 5jj' (Sai-i + (5„'+i) 

r (u)(i'j') = 7<W%'(<W Tn lj + $iL x /™>L m j)- (4) 

In the steady state (d(q n qi)/dt) = 0. From this and using 
Eq. we get the matrix equation |I8| 



ab + ba T = d = 




e 



(5) 



where b is the correlation matrix with elements b„; = 

(2T L S il +2T R S iL J. Invert- 
ing this equation one obtains b and thus all the moments 
which includes the local temperatures Tij = (p^/m^) 
and currents = {xi-ijpij/niij). The dimension of the 
matrix which has to be inverted is iV(2iV+f ) xiV(2iV+f ) 
and using the fact that it is a sparse matrix we have been 
able to numerically 0| obtain b for system sizes up to 
L x = Ly = L = 8. 

The molecular dynamics simulations were performed 
using a velocity- verlet scheme [2(|. We chose a step- 
size of At — 0.005 and averaged over I0 8 time steps (for 
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FIG. 2: Plot of disorder-averaged-current versus system size 
for the two different heat baths and for the case of correlated 
disorder with white noise. For the full disorder cases, the solid 
lines are fits to the last three points and have slopes 0.59 and 
0.51. For the case of correlated disorder, the slope from exact 
numerics (and also simulations) is compared to 1.5 which is 
what one expects analytically. 

L = 256 we took At = 0.02 and 10 7 - 5 x 10 7 time steps). 
The temperatures at the two ends were set to Tl = 0.5 
and Tr — 2.0. In Fig.QJ we plot at every site, as ob- 
tained from the simulations and from the exact solution, 
for a particular realization of disorder. The agreement is 
clearly very good. We also find that that current fluc- 
tuations decay faster than fluctations of the local tem- 
perature. This is because, in the harmonic model, decay 
of fluctuations takes place only through coupling to the 
reservoirs and this is weak for localized modes that con- 
tribute to the temperature. This means that equilibrated 
values for the current can be obtained in smaller simula- 
tion runs. 

Simulations were performed for sizes L x = L y = L = 
4, 8, 16, 256 and in Fig. |2| we plot the system size- 
dependence of the current, averaged over 100 samples 
(9 samples for L = 256). Error-bars shown are those 
calculated from the disorder averaging, the thermal ones 
being much smaller. For larger system sizes we find that 
we need to average over a smaller number of realizations 
since the rms spread in the current decreases rapidly. 
From our data we estimate a = 0.59 ± 0.01. 

We briefly note that Eq. \5\ can be solved exactly for 
the ordered case, using methods similar to those in (lj| . 
The current is independent of system size and given by 

j _ (T L T R ) I d 0^ 
4tt7 J 

where M<l) = 1 + \{l~ 2 + Q ~ \[^{l~ 2 + Q + il~ 2 + 
Iq) 2 ] 1 / 2 and I = 2[1— cos(g)]. The temperature in the bulk 
of the system takes the constant value T = (Tl + Tr)/2. 

Simulations with correlated noise: In this case 
the simulations were done by using a slightly modified 
version of the velocity-verlet algorithm with a step-size 



At = 0.001 and averaging over 10 s time steps. The ac- 
curacy of the algorithm was tested in one-dimensions 
where exact numerical results are available [f|. Simu- 
lations were performed for sizes L = 4,8, 16, 128 with 
disorder averages over 100 samples (22 for L = 128). The 
results are plotted in Fig.[21and we estimate the exponent 
a = 0.51 ± 0.01 in this case, which is somewhat different 
from the slope of a « 0.59 for the case of uncorrelated 
noise. It is possible that the small difference is a finite 
size effect and for larger system sizes we might see the 
same exponent. The error bars given are statistical er- 
rors while those from finite-size effects are more difficult 
to estimate. The next example throws some light on this 
aspect. 

Correlated disorder: Finally we consider a special 
case of correlated disorder (with white noise baths) which 
was discussed in 0. This case is analytically tractable 
and gives us some insights on possible finite size effects. 
In this model, in a given column, say the ith, all parti- 
cles have the same mass rrii . This case can be reduced to 
an effective one-dimensional problem Q. Using the fact 
that there is order in the transverse direction, we trans- 
form to new variables using an orthogonal basis ipj(q) 
which satisfies the equation 2tpj (q) — ipj- 1 (q) — tpj+i (q) = 
X q ipj(q). We choose the tpjil) to be real and find that 
X q = 2(1— cos (q)) with q — 2sir/L y where s = 1, 2, L y . 
The new variables Xi(q) — J2j x i,j' l Pj( < l) satisfy the fol- 
lowing equations of motion: 

m\Xi(q) = -pt(q)xi+x 2 -jx 1 +r/ L (q) 
niL x x Lx (q) = -fi(q)x Lx + x Lx -\ - 7±l* + V R (q) 
miXi(q) = -[i(q)xi + Xi-i + x i+ i (6) 

(for 1 < i < L x ), where /z(g) = 2 + l q . The trans- 
formed noise variables n L (q,t) = J2j r lf(t)' l l J j('l) satisfy 
(v L (q,t)v L (q\t')) = 2T Ll 5(t - t')5 qq/ and similarly for 
ij R (q,t). Thus for every q we have an equation identical 
to that of a one-dimensional disordered chain with an ad- 
ditional on-site potential V = l q x 2 /2. The heat current 
in terms of the transformed variables is: 

•/=^E«-^)+^))±i(<z)> (7) 

« q 

Fourier transforming Eq. using: Xi(q,0) = 
J dtx l (q,t)e- i0t , fj L > R (q,0) = j dtn L > R (q,t)e-' t0t we get 
the set of equations: 

[n(q) — mi0 2 + i70]ii - i 2 = V L 
—Xi-i + \p(q) - mi0 2 ]xi - x i+ i = 
-x Lm -i + [n(q) - m L;c 2 + i70]x ix = fj R 

which in matrix notation can be written as 

y(q,0)x(q,0) = f](q,0) where J, fj are the vectors 

(x 1 (q,0),...,x Lx (q,0)) T and (i) L (q, 0), 0, 0, fj R (q, 0)) T . 

The matrix y — fc$ - 2 M + i0T where 
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M n i = 5 n im n ; = fi{q)S n i - (5 ni _i - S n i +1 and 

r„; = 7^rii(^ni/ TO i + ^nL x /m Li ). After some manipula- 
tions the current in Eq. [3 simplifies to give: 

j = l2{Tl KLy TR) Y. fj^ 2 \y;lM\ 2 (8) 

The inverse element is given by 0: tyil (0, g) | 2 = 
\Det[y}\- 2 where Det[y] = D hLx + i-y0(D 2 , L , + 
£>i,l x -i) — 7 2 2 £>2,l x -i and D^i is defined to be the de- 
terminant of the submatrix of fc$ — 2 M beginning with 
the ith row and column and ending with the i'th row 
and column. These matrix elements can be expressed in 
terms of products of random matrices 0. Using these 
results one can very efficiently compute the integral in 
Eq. |H1 and obtain J accurately for quite large system 
sizes (L = 512). In Fig.[2we show the system size depen- 
dence of the current as obtained from the exact numerical 
method and also from simulations. They agree very well 
and give a rs 1.5. This value can be understood ana- 
lytically by noting that the leading contribution to the 
current in Eq. [S] comes from the q — > term (finite q 
modes decay exponentially with system size) and this is 
identical to a pure ID chain for which a = 3/2. 

The fact that the simulation results agree extremely 
well with the exact numerical results (for sizes up to L — 
128) proves the accuracy of our simulations. Further we 
see that for the correlated disorder case the asymptotic 
result for the exponent can already be seen at around 
L = 512. This gives us confidence that for the case of 
the fully disordered lattice we might already be close to 
the asymptotic value. This is also supported by the fact 
that the change in slope of the J-vesus-L curves in Fig. [21 
over the system sizes studied is very small. 

Conclusions: We have performed extensive simula- 
tions of heat conduction in a mass-disordered harmonic 
solid in 2D which give exponents a ~ 0.59 for white 
noise heat baths and a ~ 0.51 for correlated noise baths. 
A system with correlated disorder gives, somewhat sur- 
prisingly, a larger exponent a = 3/2. The combination of 
simulations and exact numerics gives us confidence on the 
accuracy of our results and also additional insight. Some 
interesting open problems are the exact determination 
of the exponent a in 2D, for any heat bath model, and 
an analytical understanding of dependence of a on bath 



properties. 
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